物种多度的非约束排序中被动添加环境变量概述

前篇先后简介了主成分分析（PCA）、对应分析（CA）、主坐标分析（PCoA）以及非度量多维尺度分析（NMDS）。这些排序方法均属于非约束排序，只涉及一个数据矩阵，并在低维空间中尽可能呈现原始的数据结构。计算过程中不存在解释变量。

在基于距离的排序方法PCoA和NMDS中，前篇以物种多度数据的排序为例，提到了物种变量信息会在计算样方距离的过程中丢失。并提到若期望补充物种得分，物种变量可以直接与排序轴进行相关分析，或者通过多度加权平均方式排序轴建立关联而投影到排序空间中。类似这个方式，环境变量（可假定为解释变量）也可被动投影到排序空间。

在群落数据的排序分析中，若我们想在排序图中尽量展示最大的物种矩阵的变化量（variation）的同时，又想展示出环境因子（即解释变量）；或者当我们在约束排序中发现环境因子的解释程度偏低，约束排序模型不是很适用时，不妨考虑试一下在非约束排序模型中被动添加环境变量的方式。

在非约束排序中，非约束轴的生态学意义可以通过将这些轴上的样方得分（sample scores）与外部变量（supplementary variables，在群落数据分析中，通常为测量或估计的环境变量）建立联系来解释。这种联系可以通过计算环境变量与前两个或几个主轴间的相关性（常用Pearson相关系数），或者通过使用（加权）多元回归进行拟合来实现。相较之下，相关系数的方法更为直观且易于理解，但仅适用于线性排序的方法（PCA等）；（加权）多元回归虽不直观，但更加通用，适用于包含线性（PCA等）、单峰（CA等）等排序方法在内的多种模型。二者差异的主要原因在于，在线性模型的排序方法中，所有样方是等权重代入计算的；而在单峰模型中，样方权重则与该样方中物种丰度总和成比例。因此，当计算外部环境变量与排序轴的关联程度时，必须充分考虑这一点，且在计算中需要将样方权重包含在内，这也是在非线性模型的排序中使用多元回归的原因。通过这些方法，能够将环境变量被动地投影在已存在的非约束排序空间中。对于环境变量与对应排序轴的关系强度，若基于相关系数的方法，则由相关系数（r）决定；若为多元回归的方法，则由回归决定系数（R²）决定；同时检验相关系数或回归决定系数的显著性也是不可或缺的。

需要注意，这种在物种多度的非约束排序模型中被动投影环境变量，以对排序轴作出解释的方法，与约束排序相比是有本质区别的。由于外部环境变量只能在非约束排序空间已经生成后，通过被动的方式投影到排序空间中，因此非约束排序本身的运算过程没有受到外部环境变量的影响，其主要表达的仍然是数据矩阵中对象（样方）或变量（物种）间的关系。在将环境变量投影至包含样方-物种信息的非约束排序图中时，此时的排序轴代表一个虚拟的环境梯度，被动加入环境因子的目的在于考察这个虚拟的环境梯度跟当前所测的环境变量之间的关系。而约束排序在计算时则明确地探索了响应变量矩阵和解释变量矩阵的关系，例如在生态群落的RDA分析中，所得样方和物种排序已经是结合了环境因子的线性组合（即是先跟环境因子进行回归的拟合值的排序）。

以下将分别简述这两种通过“被动投影”的方式，在物种多度的非约束排序中加入外部环境变量的方法。

环境变量与非约束轴的相关性

以tb-PCA为例简要描述其步骤：

（1）基于物种多度数据计算排序轴上的样方得分，即首先需执行tb-PCA分析获得样方排序坐标。（图示a-b-c）

（2）在tb-PCA结果中选择主要的排序轴，通常为前2-3个（一般不考虑很多的轴，毕竟排在后面的轴的特征根很低，不具代表性），并计算每个选择的轴的样方得分与外部环境变量之间的相关系数（Pearson相关），同时还需检验相关系数的显著性仅保留显著的相关关系。（图示c-e-f-g）

（3）将每个环境变量与选定排序轴的样方得分的相关系数，呈现为向量的形式，向量的角度和方向取决于相关系数的值和符号。例如在常见的二维坐标中，由坐标[0,0]起始，指向[r₁,r₂]。（图示g-h）

（4）将上述相关性（向量）映射到最初的排序空间中，即可在原初的物种非约束排序图中添加了环境变量信息。这样即使是PCA，也同时涵括了样方-物种-环境变量3组信息，形似RDA的效果（但和RDA有本质区别）。（图示d-h-i）

环境变量与非约束轴的（加权）多元回归

在多元回归中，将所选定的排序轴的样方得分作为解释变量（自变量），各环境变量作为响应变量（因变量），则对于每个环境变量（env_i），遵循以下等式：

其中，b₁、b₂... b_n等是回归系数，score₁、score₂... score_n等是所选定的第n个排序轴上的样方得分。由此得到的回归系数，经过标准化并乘以回归决定系数（R²）后等价于上文部分所描述的相关系数，即可以向量的形式投影在排序图中。当然，结合上文的内容我们可知这种“等价”仅适用于线性模型的排序方法（PCA等），因为在单峰模型的排序方法（CA等）中，还需考虑权重信息（非线性模型中，样方权重与样方中物种丰度有关）。对于回归R²显著性的判别，可以通过蒙特卡洛置换检验（Monte Carlo permutation test）来实现，以验证各环境变量在排序轴上的拟合优度是否合理，最终只保留显著的回归关系并将相应的环境变量映射在排序图中。

更多细节提要。本篇在一开始提到，在单峰模型的排序方法中使用多元回归添加环境变量时，需要充分考虑样方权重，样方权重与该样方中物种丰度总和成比例。在多元回归中，权重用于对所有变量中心化处理，此时计算的b₀（截距）为零，因为响应变量（环境变量）和解释变量（即所选定的排序轴的样方得分）都是经过中心化处理后的（具有零均值，但不是标准化）；回归系数b₁、b₂... b_n等需要再经过标准化处理，以使它们的平方和为1。标准化后对应的值c₁、c₂... c_n等，称为余弦，其计算公式：

其中k为比例因子（scaling factor）：

对于此时所得余弦c₁、c₂... c_n等，再乘以回归决定系数（R²），则它们将等同于前述Pearson相关系数中的r₁、r₂…r_n（同前所述，仅局限于线性回归且不考虑权重的情形中；r_n即为该环境变量在所选定的第n个排序轴中，与该轴的样方得分的相关系数）。

例如下表为两个环境变量的结果示例。PC1和PC2中的值不是相关系数，而是余弦（见上文描述，它们的平方和等于1）。R²为每个环境变量的（加权）多元回归的决定系数，该值越高，代表所对应的环境变量越具解释性。Pr (> r) 为蒙特卡洛置换检验的结果，用于判断回归变量的显著性，以帮助我们选择有效的环境变量加以分析。

相关系数或回归系数的显著性检验

无论计算相关系数或使用多元回归，关联程度的显著性检验均是不可缺少的。这里使用一个示例简要说明，如下图所示，使用多元回归在PCA排序图中被动地加入环境变量。其中，左图仅为一组随机变量，右图为真实的环境数据。

对于左侧的图来讲，由于环境变量仅仅为随机生成的一些变量，而非真实的数据，因此理论上并不应当与排序轴产生显著的拟合R²。我们可以直接在结果中观测到，所有随机变量的回归决定系数R²的值均很低，与此相对应，几乎不存在显著的p值（默认显著性p < 0.05）。但是若我们忽略回归的p值，则将缺乏选择的依据，这些低拟合的环境变量加入至排序图中参与解释的合理性就无法知晓。且在多元回归中，由于变量向量的长度取决于R²，即具有较高R²的环境变量会在排序图中产生相对较长的投影向量；即便所有变量均为随机产生，但它们仍然会具有非零R²（尽管比较低），并且可能会出现某些相对高一些的R²使得其对应的随机变量的向量长度略长。此时仅在排序图中观测向量投影时，也将会错误地认为某些投影长度较长环境变量（随机变量）具有很好的解释程度。

对于右侧的图来讲，由于加入的为真实的环境变量数据，因此其与排序轴的回归决定系数R²明显高于左图中使用随机变量的结果。具有高拟合R²的结果对应了更为显著的p值，这是毋庸置疑的，低拟合的环境变量可通过p值加以判断并剔除。

图注：使用多元回归，将9个随机变量（左图）和实际环境变量（右图）被动投影至PCA排序图中（选择前两个主轴），并执行蒙特卡洛置换检验回归的显著性。若忽略检验的显著性，则我们可能会试图解释排序图中所有的变量，即便它们可能仅仅来自一组随机数所组成；而实际上，对于随机变量来讲，通过回归不应当与排序轴间产生可信（拟合度较好）的R²；即便对于真实的环境变量数据，回归结果也并不完全都显著，即并非所有环境变量都能合理地参与解释。

因此，检验回归拟合优度的显著性显得尤为重要，以帮助我们识别有效（R²绝对值足够高）的结果，进而选择合适的环境变量投影至物种多度的非约束排序空间中参与解释。

并且，以上提及的方法（无论相关性或多元回归），主要用于反映环境变量和选定非约束轴之间的关联信息，在检验这种关系的显著性时，检验的零假设是“环境变量和选定非约束轴之间没有关联”，前提条件要求二者相互独立。

P值校正

由于在置换检验中，最小p值取决于置换次数，当存在较多的外部环境变量有待被考虑时，可能需要增加置换的次数以降低最小p值，这是很必要的。例如，当置换次数设置为99时，则置换检验可以达到的最小p值“p_min = 1 / (99 + 1) = 0.01”，若我们此时以p < 0.01作为拒绝零假设的依据，那么显然是不可能的；适当增大置换次数，例如为499次时换检验可以达到的最小p值“p_min = 1 / (499 + 1) = 0.002”，此时若以p < 0.01作为拒绝零假设的依据是可行的。

因此在使用置换检验时，置换次数需随变量数量的增加而增大，过少了肯定不合适，但同样地，次数过多也并不是一个合理的方式。为了使结果更准确，我们通常会执行更多次数的检验过程，在增加置换次数尽可能排除不可信（不显著）的结果的同时，随着检验的次数的增多，也会观测到更多的显著性结果，即使真实情况中其零假设为真（本身不显著，但随着检验次数的增多被判断为显著，即拒绝了零假设，产生I类错误，可以理解为假阳性）。

这里通过一个示例加以说明。随机生成两组正态分布的随机变量，对它们计算回归并执行参数F检验（parametric F-test），结果如下图所示。由于变量是随机生成的，理论上不会产生显著的回归结果。但在重复100后，却产生了显著的结果。显著结果的比例取决于所选取的显著性阈值，例如，将p值5％（alpha = 0.05）判断为显著的结果，那么即使变量是随机的，大约仍会存在5％的检验结果是显著的（即这5%的结果拒绝了零假设“变量间不存在关联”，由于零假设为真，产生I类错误）。换句话说，至少有一个检验结果在p < alpha时显著的概率可以表示为“1 - (1 - m)^alpha”，又称I类错误率，若我们尝试解释未经任何校正的多重检验的结果时可能会犯I类错误。

一个常用方法是在检验后加入p值校正过程，可以降低犯I类错误的方法，例如Bonferroni校正、FDR校正等，但也可能与之并存随检验次数增加而功效降低的情形，II类错误随之增加（零假设“变量间不存在关联”为假，即变量间实际上是存在某种联系的，但我们接受了错误的零假设，导致II类错误，可以理解为假阴性）。权衡之下，我们还是倾向于尽可能降低I类错误的产生，即尽可能保证所得结果都是真阳性的。

图注：生成两组正态分布的随机变量并使用parametric F-test检验回归的显著性，重复100次，每次产生新的随机变量。图中将显著（p < 0.05）的回归标记为红色，在总共100次的分析中，4次结果显著（≈5％出现次数）。

在上文阐述“相关系数或回归系数的显著性检验”的示例中，对于9个随机变量在PCA排序轴中的多元回归，执行蒙特卡洛置换检验后初步排除了绝大多数回归结果不显著的数据。但是我们需要注意，此时仍存在一个随机变量与PCA前两排序轴中的样方得分具有显著的关联，事实真的是这样吗（对于随机变量而言，“二者无关联”的零假设更贴切真实情况，此处拒绝零假设的判断是否有误）？此外，对于真实的环境变量数据的多元回归结果，是否同样存在着假阳性呢（事实上某环境变量与排序轴并无显著关联，但检验时产生了误判，拒绝了“二者无关联”的零假设，由于该零假设为真，故产生了I类错误）？

针对这一问题，引入Bonferroni校正原始所得蒙特卡洛置换检验的p值，如下图所示。此时我们可以看到，所有随机变量均不可通过PCA排序轴参与解释（包括此前未通过Bonferroni校正的结果中判断为显著的随机变量“rand 1”），这正是我们所期望的。同时，对于真实环境变量数据的多元回归结果来讲，经过Bonferroni校正后环境变量“COVERE3”不再显著（此前为显著）。

图注：使用多元回归，将9个随机变量（左图）和实际环境变量（右图）被动投影至PCA排序图中（选择前两个主轴），执行蒙特卡洛置换检验回归的显著性，并通过Bonferroni校正p值。与上文未经过p值校正的结果相比，左图随机变量全部被舍弃，右图中对环境变量做了更谨慎的取舍（暂且忽略II类错误的情形）。

p值校正后，相较于原始p值的数值会升高，这个过程会影响到置换检验的最小p值。例如当置换次数设置为99时，则置换检验可以达到的最小p值“p_min = 1 / (99 + 1) = 0.01”，理论上此时是可以将p < 0.05作为拒绝零假设的依据。对于10个补充变量，执行99次置换检验且引入p值校正过程，多重检验的p值校正方法均为Bonferroni（p值×置换次数），那么能达到的最小校正p值将是“0.005×10 = 0.05”，此时我们将无法再根据p < 0.05拒绝零假设。因此，在执行p值校正前，还需记得适当增加置换检验的次数以合理降低最小p值达到判断显著性的目的。

参考资料

主要翻译自：https://www.davidzeleny.net/anadat-r/doku.php/en:suppl_vars

友情链接